#include <math.h>
#include <string.h>

#include "eno.h"
#include "fluids.h"

#include "eno.xtn"
#include "fluids.xtn"

#define jet_wave_number 0.30
#define wave_amplitude 1.0
#define wind_speed_nomore 0.0

#define legendre_cnt 2

extern void init_stokes_data();
extern void init_wavedata();
extern double get_wavedata();
extern double compute_heaviside_general();

double density_below,density_above,density_middle,
       laplace_above,laplace_below,laplace_middle;

double hc_save,max_wave;
logical hc_init;

#define dr 0.01

double try_h(h0,rc,e0,rn)
double h0,rc,e0,rn;
{
double h1,h2,temp,temp1,rr,hprime;

  h1=h0+dr*(-rc/sqrt(1-sqr(rc)));
  h2=h0;
  rr=rc+dr;
  while (rr<rn)
   {
    hprime=(h1-h0)/dr;
    temp=sqrt(1+sqr(hprime));
    temp1=-hprime*sqr(temp)/rr-e0*h1*temp*sqr(temp)/4.0;
    h2=sqr(dr)*temp1-h0+2.0*h1;
    h0=h1;
    h1=h2;
    rr+=dr;
/*printf("rr,h2%10.5f%10.5f\n",rr,h2);*/
   }
  return h2;
}

#define rlast 2.0
 
double init_meniscus(rc,e0)
double rc,e0;
{
double h,ha,hb,a,b,c,hc;

  if (hc_init)
   return hc_save;

  a=0.0;
  b=0.5;
  c=b;
  ha=try_h(a,rc,e0,rlast);
  hb=try_h(b,rc,e0,rlast);
  hc=hb;
  while (fabs(hc)>0.00001)
   {
/*
printf("a,b,c,ha,hb,hc%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n",
  a,b,c,ha,hb,hc);
*/

    c=b-hb*fabs(b-a)/fabs(hb-ha);;
    hc=try_h(c,rc,e0,rlast);
    if (hc>0.0)
     { b=c;hb=hc; }
    else
     { a=c;ha=hc; }
/*
printf("after a,b,c,ha,hb,hc%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n",
  a,b,c,ha,hb,hc);
*/
   }

  hc_save=c;
  hc_init=TRUE;
  return c;
}

double compute_meniscus(hc,rc,xx,e0)
double hc,rc,xx,e0;
{
  if (xx>=1.0)
   return 0.0;
  else
   return try_h(hc,rc,e0,xx);
}

double compute_stream(xx,yy)
double xx,yy;
{
  if (xx<0.0)
   return 0.0;
  else if (xx<2.0)
   return 10.0*xx;
  else return 20.0;
}   /* compute_stream */

double compute_eta(xx)
double xx;
{
double temp;
  temp=sin(10*pi*xx/10.0)+0.5*sin(2*pi*xx/10.0);
  temp=get_wavedata(xx); 
  return temp/wave_amplitude;
}

double max_a()
{
int ix;
double tmax,tmin;

  tmax=tmin=compute_eta(5.0);
  for_all_x
   if (pointx_ok)
    {
     tmax=max(tmax,compute_eta(x[ix]));
     tmin=min(tmin,compute_eta(x[ix]));
    }
  return (tmax-tmin);
}
    


/* compute initial values (depending on data read from eno.dat)
    input xx,yy
    output u,v,d1,d2,d3
*/
void u0(u,v,d1,d2,d3,xx,yy)
double *u,*v,*d1,*d2,*d3,xx,yy;
{
double interface,bubble_y1,bubble_y2,bubble_d1,bubble_d2,
       rr,closest_y,bubble_rr,cosine_value,epsilon,e0,rc,phic,hc,
       zc,xc,slope,a,b,L,g,Uend,Tend,Ustart,A,B,r1,legend0,legend1,
       legend2,distance_to_corner,distance_to_circle,distance_to_surface,
       m_height,y_cutoff,xctr,y1,y2;
int i,cnt;

  *u=*v=*d1=*d2=*d3=0.0;

  if (fluid_select==91)
   {
    *u=*v=0.0;
    *d1=sqrt(sqr(xx-5.0)+sqr(yy-5.0))-1.0;
   }
  else if (dambreak())
   {
    *u=*v=0.0;
    xctr=5.0;y1=1.25;y2=3.75;
    if ((yy<=y1)||((yy<=y2)&&(xx<=xctr)))
     {
      if (xx>xctr)
       *d1=yy-y1;
      else if (yy>y1)
       *d1=max(yy-y2,xx-xctr);
      else
       *d1=max(yy-y2,-sqrt(sqr(xx-xctr)+sqr(yy-y1)));
     }
    else if (xx<xctr)
     *d1=yy-y2;
    else if (yy<y2)
     *d1=min(yy-y1,xx-xctr);
    else
     *d1=min(yy-y1,sqrt(sqr(xx-xctr)+sqr(yy-y2)));
   }
  else if (marcus_flag)
   {
    interface=2.5+compute_eta(xx);
    *d1=yy-interface;
    bubble_y2=(0.5+compute_heaviside_general(*d1-max_wave/2.0,max_wave/h));
    bubble_y1=wind_speed_nomore*bubble_y2;
    hc=(1.0-bubble_y2)*(compute_eta(xx+h/2)-compute_eta(xx-h/2))/h;
    rc=sqrt(1.0+sqr(hc));
    *u=bubble_y1/rc;
    *v=bubble_y1*hc/rc;
   }
  else if (fluid_select==80)
   {   /* Zalesaks problem */
    A=pi/314.0;
    *u=-A*(yy-50.0);
    *v=A*(xx-50.0);

    if (zalesakcmp)
     { *u=*v=0.0; }

    *d1=sqrt(sqr(xx-50.0)+sqr(yy-75.0))-15.0;
    if ((xx>=47.5)&&(xx<=52.5))
     {
      if (yy<=60.0)
       *d1=( (xx<50.0) ? sqrt(sqr(yy-60.0)+sqr(xx-47.5)) :
          sqrt(sqr(yy-60.0)+sqr(xx-52.5)) );
      else if (yy<=85.0)
       *d1=min( ( (xx<50.0) ? (xx-47.5) : (52.5-xx) ) , 85.0-yy );
      else if ((yy<=90.0)&&(*d1<=0.0))
       *d1=max(*d1,85.0-yy);
     }
    else if ((*d1<0.0)&&(xx<47.5))
     {
      if (yy<=85.0)
       *d1=max(*d1,(xx-47.5));
      else
       *d1=max(*d1,-sqrt(sqr(xx-47.5)+sqr(yy-85.0)) );
     }
    else if ((*d1<0.0)&&(xx>52.5))
     {
      if (yy<=85.0)
       *d1=max(*d1,(52.5-xx));
      else
       *d1=max(*d1,-sqrt(sqr(xx-52.5)+sqr(yy-85.0)) );
     }
   }    /* zalezak's problem */ 
  else if ((fluid_select>80)&&(fluid_select<=83))
   {
    *u=*v=0.0;  /* see solid_body_flow() for details */
    *d1=sqrt(sqr(xx-50.0)+sqr(yy-75.0))-15.0;
    if (fluid_select>81)   /* rescale */
     *d1=sqrt(sqr(xx-0.5)+sqr(yy-0.75))-0.15;
   }
  else if (rayleigh_taylor() )
    {
     *u=*v=0.0;
     interface=0.05*cos(2*pi*xx);
     *d1=yy-interface;
    }    /* incompressible rayleigh taylor problem */
  else if ((fluid_select==47)||(fluid_select==147))
   {
    *u=*v=0.0;

    *d1=sqrt(sqr(xx)+sqr(yy-1.1))-1.0;
    *d3=yy;
    *d2=-min(*d1,*d3);
    if (fluid_select==47)   /* not 3 fluid version */
     *d1=-*d2;

/*
    *u=(compute_stream(xx+h/2,yy+h/2)+compute_stream(xx-h/2,yy+h/2)-
        compute_stream(xx+h/2,yy-h/2)-compute_stream(xx-h/2,yy-h/2))/(2*h);
    *v=-(compute_stream(xx+h/2,yy+h/2)-compute_stream(xx-h/2,yy+h/2)+
         compute_stream(xx+h/2,yy-h/2)-compute_stream(xx-h/2,yy-h/2))/(2*h);
*/

/*  printf("A,B,Tend,Ustart%10.5f%10.5f%10.5f%10.5f\n",A,B,Tend,Ustart);*/
   }
  else if (drop_jet() )
   {
    *u=*v=0.0;
    epsilon=0.05;
    interface=1.0+epsilon*cos(jet_wave_number*yy);
    *d1=xx-interface;
   }
  else if (stratified_flow() )
   {
    epsilon=0.1;
    bubble_y2=-5.5+epsilon*cos(pi*xx);  /* water/oil interface */
    bubble_y1=-4.0+epsilon*cos(pi*xx);  /* oil/air interface */
    bubble_d1=bubble_y1-yy;
    bubble_d2=yy-bubble_y2;
    *u=*v=0.0;
    *d1=bubble_d1;
    *d2=bubble_d2;
    *d3=-min(*d1,*d2);
   }
  else if (drop_merge() )
    {
     bubble_y1 = -4.5;
     bubble_y2 = -0.5;
     bubble_d1 = sqrt(sqr(xx)+sqr(yy-bubble_y1))-1.0;
     bubble_d2 = sqrt(sqr(xx)+sqr(yy-bubble_y2))-1.0;
     *d1=min(bubble_d1,bubble_d2);
     if (three_phase())
      {
       if (feed_back() )
        {
         *d1=bubble_d2;
         *d3=bubble_d1;
         *d2=-min(bubble_d1,bubble_d2);
        }
       else
        {
         *d2=bubble_d1;
         *d3=bubble_d2;
        }
      }   /* three_phase */
     *u=*v=0.0;
    }
  else if (drop_periodic() )
    {
     /* P2=1/2(3x^2-1) for axisym and 1/2(2x^2-1) otherwise.
        we assume axisymmetric, gravity is the amplitude. */
     rr=yy/sqrt(sqr(xx)+sqr(yy));   /* cos(theta) - inverted in LM */

     legend0=1.0;
     legend1=rr;

     for (i=2;i<=legendre_cnt;i++)
      {
       legend2=( (2.0*i-1)*rr*legend1-(i-1.0)*legend0 )/i;
       legend0=legend1;
       legend1=legend2;
      }

     *d1=sqrt(sqr(xx)+sqr(yy))-1.0-fabs(gravity)*legend2;
     *u=*v=0.0;
    }
  else if (bubble_merge() )
    {
     bubble_y1 = 0.5;
     bubble_y2 = -2.80;
     bubble_d1 = sqrt(sqr(xx)+sqr(yy-bubble_y1))-1.5;
     bubble_d2 = sqrt(sqr(xx)+sqr(yy-bubble_y2))-1.0;
     *d1=min(bubble_d1,bubble_d2);
     *u=*v=0.0;
    }    /* incompressible bubble merge problem */
  else if ( ((oil_merge()) && (fluid_select==50)) ||
            (fluid_select==51) )
    {
     bubble_y1 = 0.5;   /* air bubble */
     bubble_y2 = -2.25; /* oil bubble */
     bubble_d1 = sqrt(sqr(xx)+sqr(yy-bubble_y1))-1.0;
     if (fabs(xx)<0.25)
      bubble_d2=fabs(yy-bubble_y2)-1.0;
     else if (xx>=0.25)
      bubble_d2 = sqrt(sqr(xx-0.25)+sqr(yy-bubble_y2))-1.0;
     else
      bubble_d2 = sqrt(sqr(xx+0.25)+sqr(yy-bubble_y2))-1.0;

     *d1=bubble_d1;
     *d3=bubble_d2;
     *d2=-min(*d1,*d2);
     if (fluid_select==51)
      *d1=min(*d1,*d3);

     *u=*v=0.0;
    }
  else if ((fluid_select==40) || (fluid_select == 43))
    {
     bubble_y1=1.0;
     if ((fluid_select==40)&&
         (fabs(density_middle-density_above)<almost_same))
      bubble_y1+=0.2;
     else if (fluid_select==40)
      bubble_y1+=0.2;
     else if (gas_burst_under)
      bubble_y1+=0.2;
 
     bubble_y2=0.0;
     *d1=sqrt(sqr(xx)+sqr(-yy-bubble_y1))-1.0;
     *d3=-yy-bubble_y2;
     *d2=-min(*d1,*d3);
     if (fluid_select==43)  /* no three-fluid here */
      *d1=-*d2;

     *u=*v=0.0;
    }
  else if (fluid_select == 44)
   {
    if (xx<0.0) xx=-xx;

    if (gravity==1.0)
     e0=4.0/surface_tension_level;
    else if (surface_tension_level==1.0)
    {
     e0=4.0*gravity;
     if (fabs(e0-4.77)<almost_same)
      e0=0.0;
    }
    else
     printf("gas burst scaling error\n");

    rc=sqrt(e0/6.0);
    phic=asin(rc);
    zc=1.0+sqrt(1.0-sqr(rc));
    hc=init_meniscus(rc,e0);
/*
printf("e0,rc,zc,hc%10.5f%10.5f%10.5f%10.5f\n",e0,rc,zc,hc);
*/

    *u=*v=0.0;
    bubble_y1=-(hc-zc+1.0);

    if (xx>rc)
     m_height=compute_meniscus(hc,rc,xx,e0);  /* positive */
    else
     m_height=hc;

    distance_to_corner=sqrt(sqr(xx-rc)+sqr(-yy+hc)); 
    distance_to_circle=sqrt(sqr(xx)+sqr(-yy-bubble_y1))-1.0;
    distance_to_surface=fabs(m_height-yy);

    if (yy>=m_height)
     {
      if (xx<=rc)
       *d1=-distance_to_corner;
      else 
       *d1=-distance_to_surface;
     }
    else
     {
      if (distance_to_circle>=0.0)
       *d1=min(distance_to_circle,distance_to_surface);
      else
       {
        if (-yy<bubble_y1)
         {
          y_cutoff=(xx-rc)*(bubble_y1+hc)/(-rc)-hc;
          if (-yy<=y_cutoff)
           *d1=-distance_to_corner;
          else
           *d1=max(distance_to_circle,-distance_to_corner);
         }
        else
         *d1=max(distance_to_circle,-distance_to_corner);
       }  /* inside of circle */
     }  /* underneath the meniscus */
    
   }  /* start pre-burst */
  else if (fluid_select==70)
   {
    *d1=sqrt(sqr(xx)+sqr(yy))-1.0;

    *u=0.0;
    *v=1.0;
   }
  else if (oil_merge() )
    {
     bubble_y1 = 0.0;    /* air bubble center   */
     bubble_y2 = -4.0;   /* oil water interface */
     *d1=sqrt(sqr(xx)+sqr(yy-bubble_y1))-1.0; 
     /* air distance */
     *d2=yy-bubble_y2;  /* water distance */
     if (*d1 <= 0.0)
       *d3=-*d1;
     else if (*d2<=0.0)
       *d3=-*d2;
     else 
       *d3=-min(*d1,*d2);

     if (fluid_select==29)
      /* oil on top of water */
      {
       bubble_d1=*d2;
       *d2=*d3;
       *d3=bubble_d1;
      }

     *u=*v=0.0;
    }
  else if (fluid_select==42)
    {
     *d1=sqrt(sqr(xx-0.88865)+sqr(yy))-.48865;
     *u=*v=0.0;
    }
  else if (fluid_select==222)
   {
    *u=sqr(xx);
    *v=0.0;
    *d1=xx;
   }
  else
    {
     *d1=sqrt(sqr(xx)+sqr(yy))-1.0;
     if (reinit_test) *d1*=1.0;
     *u=*v=0.0; 
     if (fluid_select==26)
      *v=-1.0; 
 
     if (puckett_flag)
      *d1=sqrt(sqr(xx)+sqr(yy+0.05))-0.00625;
     if (freetest_flag)
      {
       *d2=sqrt(sqr(xx-1.2)+sqr(yy+1.8))-1.0;
       L=0.5;
       bubble_d2=7*h/2;
       bubble_d1=sqrt(sqr(L)-sqr(h+L-bubble_d2));
       a=-L+bubble_d2;b=-2;
       *d1=sqrt(sqr(xx-a)+sqr(yy-b))-L;
       if (yy<b)
        {
         A=h;B=-2-bubble_d1;
         if (yy<((B-b)/(A-a))*(xx-a)+b)
          *d1=sqrt(sqr(xx-A)+sqr(yy-B));
        }
       else
        {
         A=h;B=-2+bubble_d1;
         if (yy>((B-b)/(A-a))*(xx-a)+b)
          *d1=sqrt(sqr(xx-A)+sqr(yy-B));
        }
       
       if (xx<h)
        {
         if (yy<-2-bubble_d1)
          *d1=sqrt(sqr(xx-h)+sqr(yy+2+bubble_d1));
         else if (yy>-2+bubble_d1)
          *d1=sqrt(sqr(xx-h)+sqr(yy+2-bubble_d1));
         else
          *d1=h-xx;
        }
       else if (*d1<0.0)
        *d1=max(*d1,h-xx);

       *d1+=2*h;
       *d1=min(*d1,*d2);
      }
    }    /* bubbles */

 if (colortest)
  {
   *d1=( (*d1<0.0) ? -1.0 : 1.0);
   *d2=( (*d2<0.0) ? -1.0 : 1.0);
   *d3=( (*d3<0.0) ? -1.0 : 1.0);
  }
}

/* initialize data depending on fluid_select.
   initialize h,x,y,u,v,d1,d2,d3
*/
void init_u()
{
double local_gravity,tension_temp,dim_height,dim_wavelen,periodic_domain;
int ix,iy;

 farx=fary=farxy=0.0;
 top_stream_value=bottom_stream_value=0.0;  
 /* non-zero if inflow/outflow in u */
 top_stream_value_global=bottom_stream_value_global=0.0;  
 /* non-zero if inflow/outflow in u_t */
 mx_type=REFLECT_BDRY;
 my_type=REFLECT_BDRY;
 ixmin=iymin=1;
 ixmax=mx-1;
 iymax=my-1;

 if (free_boundary_flag)
  {
   printf("doing free boundary conditions\n");
   mx_type=FREE_BDRY;
   my_type=FREE_BDRY;

   farx=fary=20.0;
   if (emad_flag)
    {
     farx=1.0;
     fary=40.0;
    }
   farxy=min(farx,fary);
  }
 else if (periodic_test()) 
  {
   mx_type=PERIODIC;
   my_type=PERIODIC;
  }

 local_gravity=0.0;
 laplace_above=laplace_below=1.0;

 hc_init=FALSE;

 if (fluid_select==91)
  {
   xmin=ymin=0.0;
   xmax=10.0;ymax=10.0;
   density_above=density_below=1.0;
   laplace_above=laplace_below=0.0;
   local_gravity=0.0;
   mx_type=my_type=REFLECT_BDRY;
  }
 else if (dambreak())
  {
   xmin=ymin=0.0;
   xmax=10.0;ymax=5.0;
   density_above=1.225;  /* air */
   laplace_above=0.0000177625;  /* air */
   density_below=999.2;  /* water */
   laplace_below=0.001137;  /* water */
   local_gravity=-9.8;
   mx_type=my_type=REFLECT_BDRY;
  }
 else if (marcus_flag)
  {
   init_wavedata();
   xmin=0.0;xmax=10.0;
   ymin=0.0;ymax=5.0;
   density_above=1.225;  /* air */
   laplace_above=1.9e-5; /* air */
   density_below=1000.0;  /* water */
   laplace_below=1.3e-3;  /* water */
   mx_type=PERIODIC;
   my_type=FREE_BDRY;
   farx=fary=farxy=20.0;
   local_gravity=-9.8;
  }
 else if ((fluid_select>=80)&&(fluid_select<=83))
  {
   xmin=-0.5;xmax=100.5;
   ymin=-0.5;ymax=100.5;
   if (fluid_select>81)   /* rescale */
    {
     xmin=ymin=-0.005;xmax=ymax=1.005;
     mx_type=my_type=PERIODIC;
    }
   density_above=density_below=1.0;
  }
 else if (rayleigh_taylor() )
   {
    if (fluid_select == 24)
     {
      mx_type = SYMMETRIC;
      xmax = 1.00;
     }
    else
     {
      mx_type = PERIODIC;
      xmax = 1.50;
     }

    my_type = REFLECT_BDRY;
    ymin=-1.0;ymax=2.0;
    xmin=0.50;
    density_above=1.0;
    density_below=7.0;
    local_gravity=1.0;

    /* below is for pucketts test problem */
    mx_type=PERIODIC;
    my_type=REFLECT_BDRY;
    xmin=0.0;xmax=1.0;
    ymin=-2.0;ymax=2.0;
    density_above=1.225;
    density_below=0.1785;
    laplace_above=laplace_below=0.0031304952;
    local_gravity=-9.80665;
   }   /* rayleigh_taylor */
 else
  {
    xmin=-2.5;xmax=2.5;
    ymin=-7.5;ymax=2.5;

    /* splash problem-drop goes down gravity is<0*/
    if ((fluid_select==147)||(fluid_select==47))
     {
      xmin=-7.4;xmax=7.4;
      ymin=-5.2;ymax=2.2;
     }
    else if ((fluid_select==44) || (fluid_select==43) ||
             (fluid_select==40))   /* gas burst problems */
     {
      xmin=-3.0;xmax=3.0;
      ymin=-4.0;ymax=8.0;
     }
    else if (stratified_flow())
     {
      xmin=-3.5;xmax=3.5;
      ymin=-7.0;ymax=0.0;
     }     
    else if ((fluid_select == 41) || (fluid_select==42))
      {
       xmin=-1.5;xmax=1.5;
       ymin=-4.0;ymax=2.0;
      }
    else if (drop_merge() )
      {
       xmin=-3.5;xmax=3.5;
       ymin=-6.0;ymax=1.0;
      }
    else if (drop_periodic() )
      {
       periodic_domain=1.5;
       if (legendre_cnt==4)
        periodic_domain=2.0;

       xmin=-periodic_domain;xmax=periodic_domain;
       ymin=-periodic_domain;ymax=periodic_domain;
      }
    else if (drop_jet() )
     {
      ymin=0.0;ymax=pi/jet_wave_number;
      xmin=-ymax/2.0;
      xmax=ymax/2.0;
     }
    else if (fluid_select==21)
      {
       ymin=-8.5;ymax=1.5;
 
       if (pressure_test)
        {
         xmin=-10.0;xmax=10.0;
         ymin=-10.0;ymax=10.0;
        } 
       else if (emad_flag)
        {
         xmin=-2.5;xmax=2.5;
         ymin=-5.0;ymax=2.0;
         if (free_boundary_flag)
          {
           xmin=-4.0;xmax=4.0;
           ymin=-9.2;ymax=2.0;
          }
        } 
       else if (puckett_flag)
        {
         xmin=-0.025;xmax=0.025;
         ymin=-0.2;ymax=0.0;
        }
       else if (terminal_flag)
        {
         xmin=-3.0;xmax=3.0;
         ymin=-10.5;ymax=1.5;
        }
       else if (axipinch)
        {
         xmin=-2.5;xmax=2.5;
         ymin=-3.50;ymax=1.50;
        }
       else if (reinit_test)
        {
         xmin=-2.0;xmax=2.0;
         ymin=-2.0;ymax=2.0;
        }
       else if (use_square_box())
        {
         xmin=-3.5;xmax=3.5;
         ymin=-5.0;ymax=2.0;
         if (reposition)
          {
           ymin=-6.0;ymax=2.0;
           xmin=-4.0;xmax=4.0;
          }
        }   /* use_square_box()==TRUE */
      }   /* fluid_select == 21 */
    else if (fluid_select==70)
     {
      xmin=ymin=-2.0;
      xmax=ymax=2.0;
     }

 
    if ((bubble_merge()) || (oil_merge()) || 
        (use_air_water()) )
      {
       density_below=density_air;
       density_above=density_water;
       laplace_below=viscosity_air;
       laplace_above=viscosity_water;
      }
    else if (fluid_select==222)
      {
       density_below=density_above=1.0;
       laplace_below=laplace_above=1.0;
      }
    else
      {
       density_below=1.0;
       density_above=weight_air/weight_helium;
      }

    local_gravity = gravity;
    laplace_middle = viscosity_oil;
    density_middle = density_oil;

    if ((fluid_select==40) || (fluid_select==43) ||
        (fluid_select==44) )
     local_gravity=-gravity;

    if ((drop_merge()) || (fluid_select == 23) ||
        (drop_periodic()) || (drop_jet()) ||
        (fluid_select==47) || (fluid_select==147))
      {
       local_gravity = -gravity;   /* was 0.0 in original test. 
                                      gravity usually 1.0 */
       density_above=density_air;
       density_below=density_water;
       laplace_above=viscosity_air;
       laplace_below=viscosity_water;
       tension_temp=tension_oil_air;
       tension_oil_air=tension_oil_water;
       tension_oil_water=tension_temp;
      }
    else if (fluid_select == 26)
      {
       local_gravity = 0.0;
       my_type = PERIODIC;
      }
   }   /* bubble problems as opposed to rayleigh-taylor */

  if (stokes_flag)
   {
    dim_wavelen=10.0;   /* 10 meters */
    dim_height=2.5;   /* 40 meters - usually */

    dim_wavelen=1.0;  /* 1 meter-- measure growth */
    dim_height=0.25;  /* no far field bc necessary here */

    if (stokes_linear_flag)
     dim_height=40.0;

    if (wind_speed>0.2)
     dim_height=40.0;   /* wind_speed acts as the slope here */

    if (stokes_breaking_flag)
     {
      dim_height=2.0*pi; 
      dim_wavelen=2.0*pi;
     }

    local_gravity=-1.0;

    density_above=density_air;
    density_below=density_water;
    laplace_above=viscosity_air;
    laplace_below=viscosity_water;
    
    xmin=-dim_wavelen/(2.0*dim_height);xmax=-xmin;
    ymin=1.0-(xmax/2.0);
    ymax=ymin+xmax;
    mx_type=PERIODIC;
    my_type=FREE_BDRY;
    fary=ymin;  /* mean nondim height is 1.0 */
    farxy=fary;
    farx=fary;   /* not used */

    if (standing_wave_flag)
     {
      dim_wavelen=2.0*pi;
      dim_height=pi;
      xmin=-1.0;
      xmax=1.0;
      ymin=-1.0;
      ymax=1.0;
      mx_type=REFLECT_BDRY;
      my_type=REFLECT_BDRY;
      if (stokes_wind_flag)
       {
        dim_wavelen=10.0;
        dim_height=500.0;
        ymin=-0.5;ymax=0.5;
        mx_type=PERIODIC;
        my_type=FREE_BDRY;
        farx=fary=farxy=dim_height-0.5;
       }
     }
    if (stokes_random)
     {
      dim_wavelen=10.0;  /* 10 meters */
      dim_height=1.25;   /* 1.25 meters - very narrow box */
      xmin=-2.0;xmax=2.0;
      ymin=0.0;ymax=1.0;
      mx_type=PERIODIC;
      my_type=FREE_BDRY;
      farx=fary=farxy=0.0;
     }
   }

  if (mx_type==PERIODIC) { ixmin=1;ixmax=mx; };
  if (my_type==PERIODIC) { iymin=1;iymax=my; };

  gravity = local_gravity;  
  /* override configured gravity except for incompressible bubble case */
  if ((use_axisymmetric == TRUE) || (use_symmetric == TRUE))
   xmin = 0.0;

  printf("xmin,xmax,ymin,ymax %10.5f %10.5f %10.5f %10.5f \n",
   xmin,xmax,ymin,ymax);
  h=(xmax-xmin)/mx;
  if (fabs(h-(ymax-ymin)/my)>almost_same)
   printf("***************deltax not same as deltay!!!! ********\n");
  else
   printf("delta x (h) is %10.6f\n",h);

  printf("initializing each grid point ...\n");

  for_all_points
   {
    x[ix]=(ix-1)*h+h/2+xmin;
    y[iy]=(iy-1)*h+h/2+ymin;
   }

  if (marcus_flag)
   max_wave=max_a();

  if (stokes_flag)
   init_stokes_data(dim_wavelen,dim_height);
  else
  for_all_points
   {
    u0(&u[ix][iy],&v[ix][iy],&d1[ix][iy],&d2[ix][iy],&d3[ix][iy],
       x[ix],y[iy]);
   }   /* looping through the grid */

}   /* init_u */

void solid_body_flow()
{
int ix,iy;
double xval,yval,xx1,yy1,A;

  for_all_points
   {
    xval=xmin+ix*h;
    xx1=xval-h/2.0;
    yval=ymin+iy*h;
    yy1=yval-h/2.0;

    s[ix][iy]=0.0;
    if (periodic_test() )
     {
      s[ix][iy]=-xval;
      u[ix][iy]=0.0;
      v[ix][iy]=1.0;
      s[ix][iy]=-xval+yval;
      u[ix][iy]=1.0;
      v[ix][iy]=1.0;
     }
    else if ( ((fluid_select==80)&&(! zalesakcmp))||(fluid_select==81))
     {
      A=pi/314.0;
      s[ix][iy]=-(A/2.0)*(sqr(xval)+sqr(yval))+50.0*A*(xval+yval);
      u[ix][iy]=A*(50.0-yy1);
      v[ix][iy]=A*(xx1-50.0);
     }
    else if (fluid_select==82)
     {
      s[ix][iy]=-sqr(sin(pi*xval))*sqr(sin(pi*yval))/pi;
      u[ix][iy]=-sin(2.0*pi*yy1)*sqr(sin(pi*xx1));
      v[ix][iy]=sin(2.0*pi*xx1)*sqr(sin(pi*yy1));
     }
    else if (fluid_select==83)
     {
      A=cos(pi*rk_time/2.0);
      s[ix][iy]=-sin(4.0*pi*(xval+0.5))*cos(4.0*pi*(yval+0.5))*A/(4.0*pi);
      u[ix][iy]=sin(4.0*pi*(xx1+0.5))*sin(4.0*pi*(yy1+0.5))*A;
      v[ix][iy]=cos(4.0*pi*(xx1+0.5))*cos(4.0*pi*(yy1+0.5))*A;
     }
   }  /* looping grid */
}   /* solid_body_flow() */
     



